/************************************************************************
*************************************************************************
*************************************************************************

simulate_va.do

Compares AR to value-added production function
		
*************************************************************************
*************************************************************************
************************************************************************/


args jobid sdset

local logon 1


clear
clear matrix
set more off


if mi("`sdset'") {
	local sdset 1
}


//Store current state of the random number generator
global curseed = c(seed)

//Load random seeds
do "do/rngseeds`sdset'.do"

//Load globals
do "constants/simglobals.do"


if $ncores >  1 {
	set seed ${seed`jobid'}
}
	

	
//Set maximum number of iterations
global prevmaxiter = c(maxiter)
set maxiter 100


if `logon' == 1 {
	cap log close

	cap mkdir "logs"
	local tmp = string(date("`=c(current_date)'", "D M Y"),"%tdCCYYNNDD")
	cap mkdir "logs/`tmp'"
	local tmp2 = clock("`=c(current_time)'","hms")
	log using "logs/`tmp'/`tmp2'simulate_va`jobid'"
}



//Targets for fraction of firms constrained
global consttargets 80 70 60 50 40 30 20 10 1




/************************************************************************
*************************************************************************
Part 1: Helper Programs
************************************************************************
************************************************************************/

//Load data generator
do "do/programs/gendatva2.do"


//Load GNR program
do "do/programs/acf2.do"



/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end





/************************************************************************
Program: truemoments
Purpose: Computes true moments (must tell it the true production function)

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop truemoments
program truemoments
	syntax name(name=name), pk(real) pl(real) [  mu(real 0) sc(real 0) ]
	
	/************************
	Compute marginal products
	***********************/
	
	/****** CES ********/
	if "`name'" == "ces" {
		//For CES they need to pass mu
		
		if mi("`mu'") {
			error "Must pass mu (elasticity of substitution) for a CES production function"
		}
		
		/* There seems to be an error here
		//Marginal product of K and L
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x''*(1-`pm') * Y^(  (1-`pm'-`mu')/(1-`pm')  ) * exp(omega + eps)^( `mu'/(1-`pm') ) / `X'^( 1 - `mu' )
		}
		*/
		
		tempvar I
		gen `I' = `pk'* K^`mu'  +   `pl' * L^`mu'
		
		foreach X of varlist K L {
			local x = lower("`X'")
			
			
			
			//This is actually the elasticity times Y/X
			gen MP_`X' = `sc' *`p`x'' * `X'^`mu' / `I'  * Y/`X'
		}
		
		di "This is CES..."
		
		
	}

	
	/****** Cobb-Douglas ********/
	
	else {
		//Marginal product is just param * Y/X
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}	
	}
	
	
	
	
	/************************
	Moments: Elasticities and Sdev of marginal products
	***********************/
	
	//Elasticity = MP_X * X/Y
	foreach X of varlist K L {
		local x = lower("`X'")
		
		//This may vary by firm
		gen elas`x'_tru 	= MP_`X' * `X'/Y 
		
		//Exclude the top and bottom 1 percent
		bysort t: cumul MP_`X', gen(MPpct)
		
		bysort t: egen sdevMP`x'_tru 	= sd(MP_`X') if MPpct > .01 & MPpct < .99
		
		//Clean up
		drop MP_`X' MPpct
	}

	
end





/************************************************************************
Program: getconstraints
Purpose: Sets a global equal to the proper values of the constraint set

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraints
program getconstraints
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist 15(1)45 {
			
			spec`spec' 2^`const'

			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end






/************************************************************************
Program: ar_dem_va
Purpose: Estimates the AR method, demeaning first
************************************************************************/
cap program drop ar_dem_va
program ar_dem_va
	preserve
	tempvar tmp
	foreach var of varlist y k l kXk lXl kXl {
		egen `tmp' = mean(`var')
		replace `var' = `var'-`tmp'
		drop `tmp'
	}
	
	#delimit;
	cap gmm (y - {bk}*k - {bl}*l - {bkk}*kXk - {bll}*lXl - {bkl}*kXl 
								- {rho}*(l.y - {bk}*l.k - {bl}*l.l - {bkk}*l.kXk - {bll}*l.lXl - {bkl}*l.kXl ) ), 
								instruments(l.(k l kXk kXl lXl) l2.k ) 
								//instruments(l.(k l m kXk-mXm) l2.m l2.k l2.kXm) 
								deriv(/bk 	= -k 	+ {rho}*l.k) 
								deriv(/bl 	= -l 	+ {rho}*l.l)
								deriv(/bkk 	= -kXk 	+ {rho}*l.kXk)
								deriv(/bll 	= -lXl 	+ {rho}*l.lXl)
								deriv(/bkl 	= -kXl 	+ {rho}*l.kXl)
								deriv(/rho 	= -( l.y - {bk}*l.k - {bl}*l.l - {bkk}*l.kXk - {bll}*l.lXl - {bkl}*l.kXl ) ) ;
	#delimit cr
	
	restore

end





/************************************************************************
*************************************************************************
Part 2: Scenarios
************************************************************************
************************************************************************/


/************************************************************************
Fixed parameters
************************************************************************/

//List of constraints (in log base 2)
//local conlist -3(2)7


//CES, not very cobb-douglas
global mu2		= .8   //Implies an elasticity of substitution equal to 5

//These are calibrated to match the ACF estimates from Gandhi et al, 2013, Industry 311, Chile
global cespk2	= .018
global cespl2	= .93
global sc2		= 1.1

////Standard deviations

//Anticipated shock
global sigeta1 	= .0925324 
global sigeta2 	= ${sigeta1} * 2

//Unanticipated shock
//global sigeps 	= 1.07
global sigeps 	= .2542193



//Number of firms
global nfirmnorm 	2613
global nfirmsmall 	1000


//Number of years
global nyearnorm 	7
global nyearshort 	4

//Persistance of productivity
global rhocal .7709392


/*********************
Production Functions
**********************/

global XKL3 "gen X_kl = (  ${cespk2}*K^${mu2} + ${cespl2}*L^${mu2} ) ^ ( ${sc2}/${mu2} )"



/*********************
Productivity processes
**********************/
global nlomproc "86.30933 - 7.081727 - 33.35803*(l.omega + 7.081727) + 4.568181*(l.omega + 7.081727)^2 - .2030311*(l.omega + 7.081727)^3 "




/************************************************************************
Specifications
************************************************************************/
local nspecs = 2

//Drop programs in memory
forvalues i = 1/`nspecs' {
	cap program drop spec`i'
}






//Baseline
program spec1
	args const
	gendatva2 "${XKL3}" `const' ${rhocal} ${nfirmnorm} ${nyearnorm} ${sigeps} ${sigeta1} 0 "" 
	truemoments ces, pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end


//Nonlinear transition function
program spec2
	args const
	gendatva2 "${XKL3}" `const' ${rhocal} ${nfirmnorm} ${nyearnorm} ${sigeps} ${sigeta1} 0 "${nlomproc}"  
	truemoments ces, pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
end









/************************************************************************
*************************************************************************
Part 3: Sweep
************************************************************************
************************************************************************/
local nreps = ${B}

cap mkdir "dta/simulations"
cap mkdir "dta/simulations/sweeps_env_va"

if ${ncores} > 1 {
	cap mkdir "dta/simulations/sweeps_env_va/`jobid'_`sdset'"
	local savefol "dta/simulations/sweeps_env_va/`jobid'_`sdset'"
	
	di "*******************"
	di "I am job `jobid'"
	di "*******************" _n _n _n
}
else {
	local savefol "dta/simulations/sweeps_env_va"
}


/*
Estimate time to complete
How many operations?
*/

local totops = `nspecs' * wordcount("${consttargets}") * `nreps'
local opscom = 0
local esttime "???"


foreach spec of numlist 1/`nspecs' {

	cap mkdir "`savefol'/spec`spec'"
	
	di _n "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"

	//Get the constraints for this specification
	getconstraints `spec'
	
	local conlist "${const`spec'}"
	local ind = 1
	foreach const of numlist `conlist' {
		di "--------------------------------------------------------------------------------------------------------"
		forvalues b=1/`nreps' {
		
			di "On Specification `spec'" _col(25) "Constraint = `const'" _col(70) "Rep = `b'" _col(85) "Est. Time: `esttime' min"
		
			quietly {
				timer on 1
		
				/********
				Generate data and compute true moments
				********/
				spec`spec' 2^`const'
				order elas* sdevMP*
				
				//Construct second-order terms
				order k l m, last
				product k l m
				
				/********
				AR Method with Translog-2
				********/
				
				ar_dem_va
				if _rc == 0 {
					gen elask_ar = _b[/bk] + 2*_b[/bkk]*k + _b[/bkl]*l
					gen elasl_ar = _b[/bl] + 2*_b[/bll]*l + _b[/bkl]*k
				}
				else {
					gen elask_ar = .
					gen elasl_ar = .
				}
				
				
				
				

				
				/********
				Ackerberg-Caves-Frazer, two-step
				********/
				if `spec' == 2 {
					cap acf2 y k l m "donl"
				}
				else {
					 cap acf2 y k l m 
				}
				if _rc == 0 {

					gen elask_acf2 = _b[k] + 2*_b[kXk]*k + _b[kXl]*l
					gen elasl_acf2 = _b[l] + 2*_b[lXl]*l + _b[kXl]*k
				}
				else {
					gen elask_acf2 = .
					gen elasl_acf2 = .
				}
				

		
				
				
				/********
				Compute estimated moments and form observation
				********/
				keep elas* sdevMP* 
				collapse _all
				
				tempfile `b'
				save ``b''
				
				
				/********
				Estimate time remaining
				********/
				timer off 1
				timer list
				local opscom 	= `opscom'+1
				local timesofar = r(t1)
				local avgtime	= `timesofar'/`opscom'
				local timerem 	= `avgtime' * (`totops'-`opscom') / 60
				local esttime 	= string(`timerem', "%9.2f" )
				
			}
			
		
		}
		
		quietly {
			use `1'
			forvalues k=2/`nreps' {
				append using ``k''
			}
			
			
			
			//if there's a negative, change it to an underscore
			gen const 	= `const'
			local fc : word `ind' of $consttargets
			gen fraccon	= `fc'
			
			gen spec = `spec'
			
			save "`savefol'/spec`spec'/const`fc'.dta", replace
			
			local ind = `ind'+1
		}
	}
	
}




//Restore default max number of iterations for GMM
set maxiter ${prevmaxiter}


if `logon' == 1 {
	cap log close
}
